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Abstract 

We present a novel approach for the inverse problem in electrical impedance tomog- 
raphy based on regularized quadratic regression. Our contribution introduces a new 
formulation for the forward model in the form of a nonlinear integral transform, that 
maps changes in the electrical properties of a domain to their respective variations in 
boundary data. Using perturbation theory the transform is approximated to yield a 
high-order misfit function which is then used to derive a regularized inverse problem. In 
particular, we consider the nonlinear problem to second-order accuracy, hence our ap- 
proximation method improves upon the local linearization of the forward mapping. The 
inverse problem is approached using Newton's iterative algorithm and results from simu- 
lated experiments are presented. With a moderate increase in computational complexity, 
the method yields superior results compared to those of regularized linear regression and 
can be implemented to address the nonlinear inverse problem. 

keywords: Impedance tomography transform, quadratic regression, Newton's metliod 



1 Introduction 

In Electrical Impedance Tomography (EIT) voltage measurements captured at the boundary 
of a conductive domain are used to estimate the spatial distribution of its electrical proper- 
ties. The technique has numerous applications in exploration geophysics ^45j, environmental 
monitoring and hydrogeophysics [s], [24], biomedical imaging [Tt], industrial process moni- 
toring (39], archaeological site assessment 34 and non-destructive testing of materials [38]. 



Owing to its many practical uses and intriguing mathematics, EIT has seen numerous theo- 
retical and computational developments, e.g. the the chapter expositions in [l], [25] and [22] , 
Among its fundamental challenges remain the nonlinearity and ill-posedness of the inverse 
problem, which inevitably compromise the spatial resolution of the reconstructed images. 
From the mathematical prospective, this inverse boundary value problem, formalized by 
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the seminal publication of Calderon presents a number of implications on the existence, 
uniqueness and numerical stability of the solution [6], [ij. Although the issues of existence 
and uniqueness can be eradicated under some mild assumptions, see for example [4l] for 
isotropic conductivity fields, the instability causes the problem to be extremely sensitive to 
inaccuracies and small errors in the data. To alleviate the ill-posedness one usually resorts in 
implementing some type of regularization strategy that stabilizes the solution \ l8j . Based on 
prior information about the unknown electrical parameters and/or the noise statistics in the 
measurements, regularization schemes are applied in order to stabilize the reconstructions. 
In the typical variational framework for example, regularization methods are often expressed 
as additive penalty terms augmenting the associated data misfit function, essentially biasing 
the solution away from features that are inconsistent with the available a priori information. 
In this sense, [i20| examines the case of Tikhonov regularization in the context of nonlinear 
system identification emphasizing the bias-variance trade-off on the solution. 

The nonlinearity inevitably increases the complexity of the problem, as the data misfit 
function has several local minima, and hence one is faced with the challenge of locating 
the solution that corresponds to the global minimum. Aside a few notable exceptions, like 
the d-bar method |26j and the factorization method [15| , algorithms that treat the nonlin- 
earity are essentially Newton-type iterative solvers, such as the often used Gauss-Newton 
(GN) method, that implement local linearization and regularization, essentially exploiting 
the Frechet differentiability of the analytic forward operator, to yield at each iteration a 



quadratic error function with respect to the unknown parameters [9], 32 . Starting from a 
feasible guess and, in some cases, an estimate of the noise level in the data, one applies a 
number of linearization-regularization cycles until a convergence is reached in the sense of 
the discrepancy principle. Analysis on the convergence rates of the GN algorithm and quasi- 
Newton variants for high-dimensional problems can be found in [s], |18| and [23| , and [14| . 
These results state that convergence is not guaranteed unless a stable Newton direction, 
descent in the usual case of minimization, is computed at each linearization point. In turn, 
this relies on the optimal tuning of regularization at each iteration, indeed a delicate and 
challenging task as the degree of ill-posedness may vary significantly |27|. To rectify this 
problem and aid convergence line search algorithms can be used, that scale optimally the 
solution increment in the descent Newton direction as indeed trust-region methods 
although more computationally complex. An additional important complication may arise 
when the typically-neglected linearization error is significantly large, invariably when the lin- 
earization point is 'not close enough' to the true solution [36]. This implies that a component 
of the linearized data should not be considered in the fitting process, since local lineariza- 
tion approximation is accurate in a rather narrow trust region, and hence to cope with the 
lack of this information at each iteration one seeks to recover a 'small' perturbation of the 
parameters. With this as background, the work in this paper focusses on the following con- 
tributions: (i) A nonlinear integral transform as a forward model that maps arbitrarily large, 
bounded changes in electrical properties to changes in boundary observations. Effectively, 
this replaces the linear approximation involving the Jacobian of the forward mapping f29|. 
The transform has a closed form and admits a numerical approximation using the finite ele- 
ment method, (ii) Exploiting the new model, a high-order misfit function is formulated for 
the inverse problem in the context of regularized regression. Numerical experiments on the 
resulting inverse problem have yield solutions with small image errors and adequate spatial 
resolution. 
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Higher-order derivatives are thus seldom used in inversion schemes since the increase in 
convergence rates may not compensate adequately for the computational effort required in 
computing the derivatives, in particular when high-dimensional discrete models with ten- 
sor parameters are concerned. Moreover, if the data misfit residual is small then the error 
contribution of the higher-order terms becomes negligibly small. The majority of inversion 
algorithms, as indeed the general theory, for nonlinear inverse problems utilize merely a first- 
order approximation of the underlying model |18| . A notable exception is the second-order 
method for nonlinear, highly ill-posed parameter identification problems in some classical 
partial differential equations, such as Helmholtz, diffusion and Sturm-Liouville jl6] . In |16j , 
the authors propose iterative predictor-corrector schemes encompassing Tikhonov regulariza- 
tion that approximate the second-order solution without solving a quadratic equation. Using 
this computationally efficient framework, they report on advantages in the final reconstruc- 
tions and significant improvements in the number of iterations required for convergence. 

As we develop our methodology we address mainly the EIT problem with complex 
isotropic admittivity and the complete electrode boundary conditions [40] . However, our 
derivations are not constrained by isotropic or complex property assumptions and thus can 
be easily shown to hold true for the similar problems of Electrical Resistance and Capaci- 
tance Tomography (ERT/ECT) with purely real coefficients in scalar or tensor field material 
properties [33j. Moreover, we show that the form of the new forward model remains un- 
changed with the governing elliptic differential equation is addressed in the context of more 
generalized boundary conditions that resemble more simplistic electrode models convention- 
ally encountered in the geophysical setting [s] , [2] . 

1.1 Notation and paper organization 

Consider a simply connected domain B c M'^, d - 2,3 with Lipschitz smooth boundary 
dB and a space depended isotropic admittivity function 7(x,cj) ■ B ^ C At an angular 
frequency u >0, the admittivity can be expressed as 

j{x,uj) - C7(x) + iuje{x.), 

where oo>Ci>a>ci>0 and 00 > C2 > e > denote the domain's electrical conductivity 
and permittivity respectively for some positive bounding constants Ci, C2, c. If there are no 
charges or sources in the interior of B and the angular frequency of the applied currents is 
small enough, then Maxwell's equations describing the electromagnetic fields in the interior 
of the domain reduce to the elliptic equation 

V • [7(x,tj)Vn(x,u;)] = 0, x e (1) 

where u denotes the scalar electric potential function. Measuring the potential at the acces- 
sible parts of the boundary of the domain through a finite number of sensors yields a set of 
observations (" that are likely to suffer from some type of noise and measurement impreci- 
sion rj. We will assume EIT systems equipped with L electrodes exciting the domain with a 
sequence of currents / = (/^, . . . , J"^), with /* = (/i, . . . , J^) all fixed at frequency 00. In such a 
case C is typically a linear combination of the electrode potentials U{P) = {Ui, . . . ,Ul), for 
i - 1, . . . ,q at the various current patterns. For an applied current pattern / we associate an 
electric potential field u in B, and an array of electrode potentials U at dB. When required 
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by the context we shall denote their dependence on admittivity and applied current as u{^) 
and u(I), or both as m(7,/); and respectively U{'y), U{I) and C/(7,/). The first and sec- 
ond partial derivatives of u with respect to 7 will be denoted by d^u and dj^u, a notation 
adopted for both continuous and discrete spatial functions, while for matrices and vectors 
the differentiation is to be considered element- wise. The position in B is specified by the 
vector X € M°', while the outward unit normal vector at the boundary is denoted n. Matrices 
and vector fields are expressed in bold capital letters while vectors and scalar fields in small 
case regular. For a matrix A, aj will denote the jth. row, Aij its (z,j)th element and A' its 
transpose. For a vector v, vi is the ith element and v is the complex conjugate. The spaces 
of real and complex numbers are given by M and C, while we use 9^{c} to express the real 
component of the complex argument c. 

The paper is organized as follows: We begin with a brief review of the the EIT model 
equations and associated preliminary concepts and then proceed to formulate the inverse 
problem commenting on existing algorithms the address the problem through local lineariza- 
tion. The next section is devoted to the derivation of the nonlinear integral transform under 
the complete electrode model and its generalization to the Poisson's equation with mixed 
boundary conditions. Further on we consider the high-order regularized regression problem 
and approximate the nonlinear system as a quadratic operator equation. Using the finite 
element we obtain a numerical approximation and subsequently implement Newton's algo- 
rithm to solve the problem. Finally, we present numerical results from simulated studies 
that demonstrate the advantages of the proposed methodology and we end the paper with 
the conclusions section. 

2 EIT model equations and preliminaries 

The complete electrode model in electrical impedance tomography is derived from Maxwell's 
time-harmonic equations at the quasi-static limit and describes the electric potential field 
in the closure of a conductive domain B with known electrical properties 7 and impressed 
boundary excitation conditions. The model has been extensively reviewed and analyzed in 
several publications, including [40] where the authors prove the existence and uniqueness of 
the solution, under some continuity assumptions on the interior admittivity. With reference 
to figure [T| assuming no charges or current sources in the interior of B, when a current 
is applied at the boundary, the electric potential u satisfies the elliptic partial differential 
equation ([T|. The applied current, inducing this field, is expressed by the Neumann boundary 
conditions 

/ ds 7(x,a;)Vu(x) -11= I^, x e Fe,, ^ = 1, . . . , L, (2) 

Jet 

7(x,a;)Vu(x) -11= 0, x e 55 \ Fe, (3) 

where Fe = UlliTe^- An accurate model of the electrodes is critical when comparing exper- 
imental measurements to synthetic model predictions. In effect, the voltage measurement 
recorded at the £th electrode with contact impedance is given by the Robin condition 

C/^ = ti(x) + 2;^7(x,cj)Vn(x) • n, x e Fg^, ^ = 1, . . . , L, (4) 

assuming that the characteristic function of the contact impedance is uniform on each elec- 
trode and > 0. The model admits a unique solution (li, U) upon enforcing the charge 
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Figure 1: The domain under consideration B with L round surface electrodes Fe^ attached 
at the boundary Fg. 



conservation principle on the applied currents and a choice of ground is made. Maintaining 
the conventional notation of ^40j, \25\ , and |6j we these constraints imply 



and / dsu = 0, 
1=1 -'^^ 



(5) 



where the applied currents should sum up to zero and the induced potential should have 
a vanishing mean on the boundary. For the so-called forward or direct problem ([l])-((5]) we 
adopt the following essential assumptions 25 , |22|. 



Assumption 1 (a) The domain B is simply connected with boundary dB at least Lipschitz 
continuous. 

(b) The electrical admittivity 7 e L°°{B) with essinf(7) > ci > 0. 

(c) The potential field u e Hl{B) = {u e H^{B) : Jq^ dsu = 0} 

(d) The applied currents I and measured voltages C, belong in the Hilbert spaces of the L, 
and respectively m dimensional complex vectors and C^, where m> L. 

We will often refer to the solution {u,U) e H^{B) ® as the direct solution, and to 
the problem Q-Q as the direct problem. Pertinent to this model is the adjoint forward 
problem [2] . Consider the direct solution under a pair drive current pattern I'^ with positive 
and negative polarity applied at electrodes Cp and e-n respectively. Moreover, let the fc'th 
boundary measurement be of the form 

Ck^Ue,-Ue,, p',n' e {!,..., L}, k^l,...,m 
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for a pair of electrodes Cp' and Cn'- In the practical setting of EIT or ERT, see for example 



the applications discussed in [17|, 13 , [24] and 43 , instead of measuring the electrode 
potentials, it is usual to measure the potential between adjacent electrodes. When captured 
systematically, this differential type of measurement yields m> L linearly independent data. 
Based on this measurement definition, the adjoint field solution {v, V) e Hl{B)®C^ satisfies 
the equations 

V- [7(x,a;)Vv(x)] = x e (6) 
7(x,a;)Vv(x) -n = x e \ T^, (7) 

ds7(x,a;)Vt;(x) -n = xeTe,, (8) 



v{x) + za{:>^,oj)\/v{x.) -n = V^, x e Te^, £ = 1, . . . ,L (9) 



where ^(x,u}) - 7(x, -w) is the conjugated admittivity, and /*" e C'^ is the adjoint current 
pattern whose ^'th entry equals to - if, if i - Cp or £ - en and zero otherwise. The 
uniqueness of the adjoint solution is subject to the constraints of the type in ([5|. 



2.1 Green's reciprocity 

In what follows, we make reference to the reciprocity principle. Originally derived from 
Maxwell's laws of electromagnetics, Maxwell's reciprocity principle has an analogue for ir- 
rotational fields known as Green's reciprocity. In the context of the impedance experiment 
it states, that if a current intensity / is applied at the boundary of a closed domain between 
two electrodes, say Pi = (ep,e„), then the potential measured at the boundary through an- 
other pair of electrodes P2 = (ep',en') will be equal to the potential measured at Pi if the 
same current is applied to P2. Impedance data acquisition instruments rely on this principle 
to avoid making redundant, i.e. linearly dependent, measurements. To see this consider a 
linear conductive medium B whose admittivity 7 has a nonzero imaginary component at the 
operating non- resonant frequency oj. Suppose we apply a time-harmonic electric current I*^ 
at the boundary of the domain through electrodes Pi, 

/'^(x,t) = J(x,a;)e*'^*, x e aP, 

where J is the current density field. From Maxwell's laws the electric and magnetic fields 
E, and H, within the domain satisfy 

VxH(x) =7(x)E(x), xeP. 

As the domain is simply connected, using E(x) = -Vu(x) and V x H(x) = J(x) reduces to 
Ohm's law 

J(x) = -7(x)Vu(x). (10) 

Let the magnitude of the applied current be equal to Iq, such that - Iq and I^^^ - -if^- 
Similarly, allow I"^ a diff'erent current pattern of unit magnitude applied through a different 
pair of boundary electrodes, say P2, inducing a new electric potential field. We denote the 
two fields as u{I'^) and u{I"^) to emphasize their dependance on the excitation currents. 
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Taking the normal component of the vector fields in (10) for multiplying with u{I^) and 
integrating over the boundary yields 

[ ds'u(/™)J(/'^)-n-- [ ds-fu{I'^)\/u{I'^)-n. 

JdB JdB 

At X € dB, let j(x) = J(x) • n be the normal component of the boundary current density 
field, then combining with conditions ^ and Q the left hand size of the equation above 
reduces to 

= tie L ds{Ue{n-za{n) 

where the last simplification follows as the supports of I'^ and I™' are disjoint. Using the 
Green's first formula, the right land side of the same equation can be developed to 

- f ds7u(/")Vu(/'^)-n = - f dx7Vu(/"^)-Vn(/'^), 

J dB J B 

and hence equating the two yields 

Working similarly for the adjoint field u(I"^) leads to 

C/e^,(/'^)-C/e„,(/'^) = - f dxjVuin-Vuil'^), 

■J B 

therefore for = 1 we have the standard form of Green's reciprocity theorem 

An alternative way to formalize this important result is via the complete electrode admit- 
tance operator A^y^z ■ C^, a complex Hermitian matrix that is the discrete equivalent to 
the Dirichlet-Neumann operator encountered at the analysis of the continuum EIT model j6| . 
For a fixed pair of 7 e L°°(B) and z e this bounded operator maps linearly the electrode 
potentials to the boundary currents inducing them, A^^zU - I. Let fJ-d, fJ'm ^ be two 
vectors of zero sum 
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and consider the current pattern I*^ = lofJ-d where "Ef^i if - 0. By the Hermitianity of A^^z 
the /c'th measurement d^k - Ue' - Ue' of the data vector C e C™' can now be expressed as 



■7,2 
-1 



fJ-d 



= f^'dUin, 



thus we arrive at the principle (11) 



2.2 The inverse problem and its linear approximation 

The inverse problem of EIT is to reconstruct the admittivity function 7 e L°°{B) given the 
operator A^^z- Invariably, this requires determining 7 given a finite set of hnearly indepen- 
dent current patters {1^,1^,...,!'^) and their respective electrode potentials {U^,U'^,. . . ,U'^). 
Typically, in EIT measurements one deals with frame(s) of (independent) data (" that arise 
as linear combinations of the U vectors. To address this ill-posed problem some prior infor- 
mation on the data noise rj and the (spatial) properties of 7 are needed. To approach this 
problem one usually considers the nonlinear operator equation 

C = f(7)+r/, (13) 

where £ ■ L°°{B) -> C". A solution to this problem can be obtained by considering the 
regularized regression problem 

7* = argmin{||C-£^(7)f + ^(7)}, (14) 

7 

where Q ■ L°° {B) M is a regularization functional. The choice of Q depends on the a 



priori knowledge on 7, and it usually takes the form of a smoothness enforcing term 37 
an LI norm allowing for sparse solutions [12] or a total variation norm that preserves large 
discontinuities in the electrical properties (7l|44|. As the forward operator was proved to 



be analytic [6], then subject to the differentiability of Q, problem (14) becomes suitable 
for gradient optimization methods 18 . Linearizing £ locally within a sphere S^p,^ - {7 ■ 



||7-7p|P < K^}, centered at an a priori guess-estimate 7p e L°°{B), yields the Taylor series 
expansion 

£h\S^„.) = £{jp) + 9^f(7p)(7 - Ip) + 0(ll7 - Ipf), (15) 

where dj£ ■ L°°{B) -> is the Frechet derivative of the forward operator and k > can 
be thought to be the Taylor series convergence radius. Truncating the series to first-order 
accuracy yields the linearized approximation of (|13|) 



C-^(7p) + 97^(7p)(7-7p) + ??, 7eS'^j„K, (16) 



which upon inserting into problem ( 14 ) leads to the regularized least-squares problem- that 
coincides with the first iteration of the regularized GN algorithm, for the optimal admittivity 
perturbation 

(57; = arg min [\\5C - d^£{^^)5^f + g{6^)\, 5C-C-£{lp). (17) 
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Given the invertibility of the Hessian ^d^£ (jp)' d^S (jp) + d^^Gijp)], implementing a GN 
algorithm for p = 0, 1, 2 . . . yields a sequence of solutions {70, 7i, 72, • • ■} that converges to 
a point in the neighborhood of 7*, subject to the level of noise in the data. Analysis 
and numerical results on the implementation of GN for the problem (17) can be found 
in many publications and textbooks on EIT, see for example [43], [I8], [27j and [23]. We 



emphasize that this popular approach, as well as its variants of Levenberg-Marquardt 18 



and quasi-Newton schemes 14 , rely fundamentally on the local linearization of the forward 
operator £, and thus yield a linear regression problem. Moreover, when Q is quadratic, 
the resulting cost-objective function to be minimized is quadratic and thus Newton- type 
methods provide for speedy analytically expressed solutions. The Noser algorithm proposed 
in [9] is a typical example of this approach, where the solution is computed after a single 
regularized GN iteration. Here we propose an alternative approach that leads to high-order 
regression problems. In this study we address explicitly the quadratic case. The starting 
point toward this direction is the nonlinear integral admittivity transform that we derive 
next. 



3 Nonlinear integral transform 
3.1 Perturbation in power 

To derive the nonlinear transform that maps changes in admittivity to those they cause on the 
observed boundary data we follow an approach of power perturbation. The method, which is 
due to Lionheart, has been developed in 37 and [36] to treat the real conductivity problem. 
Here we extend it to the complex admittivity case incorporating also the nonlinear terms 
arising in the perturbation analysis. With minimal loss of generality we restrict ourselves 
to the case of real contact impedance. If 7 and u are smooth enough, then applying the 
divergence theorem to Q for a test function e H^{B) we have 

- J^dx tpV ■ 7Vn - - J^dx ^Vu ■ Vip + J^^ ds ijj 7VM ■ n. (18) 

If Tp is set to satisfy the boundary conditions on the applied currents Q, the above becomes 

f Li L 

/ dx7Vii-VV'=X / ds {tp-^£)'yVu-n+Y,l£^i, 

where ^' e is a test vector for the electrode potentials. Plugging in the boundary condition 
on the measurements Q yields the weak form of the forward problem ^22j 



f dx7VwVV + E- f ds{tP-^e){u-Ui)^Y,Ii^i, (19) 

for all (ip,^) e H^{B) ® C^. Existence and uniqueness of the weak (variational) solution 
(u, U) e HI{B) ® has been proved in 40 . If > 0, then substituting -ip - u, - U into 
the weak form yields the power conservation law 

f dxj\Vu\^ + Y.ze f ds|7Vu-n|2 = X^^f^' (20) 
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which states that the power driven into the domain is either stored as electric potential or 
dissipated at the contact impedances of the electrodes. Consider now a complex perturbation 
7-^7 + 57, causing u ^ u + 5u in the interior, and ^ U£ + 5U£, j ^ j + Sj at the boundary. 
Recall that the normal component of the current density field at the boundary is j - jVu-n, 
under the new state of the model the volume integral in ( |20[ ) becomes 

J B J B J B 

+ / dx 5u ■ Vu + / da;7|V(5u| 

J B J B 



2 



dx 5"f\S/{u + 8u)\ 



2 



Notice that \l5u ■ Vu - Vn • Vdu hence the second and third integrals on the right sum up to 
2 dx7 *H{Vn • V5n}. For I and Z£ fixed, the surface term in (20) gbecomes 

Y^ze f ds \j + 5jf ^Y.^^ [ (bf + j dj + 5jj + \6jf), 

hence putting together the power conservation law for the new state of the model and 
subtracting (20) gives 



Y^e^Ue = J^dx jIvSu]"^ + J^dxj\/u-\/Su + J^dx6'j\\/{u + 6u)\'^ 
+ / dx 7V(5n • Vn + ^ Z£ / dsdjj + Y^t / dsj(5j 

From the weak form ( |l8| ) with ip - 6u, the second integral above simplifies as 
/ dxj\/u-\/5u - / dsSuj'Vu-n 

JB JdB 

- y~ 5uj 

L r _ 

= Y / ds{5Ui-z,5j)j 
e=i ~'^<'i 

L ^ r 

= E ^<^^^i ~ E / ds 6jj, 

thus substituting back into the previous equation gives the perturbed power conservation 
law 

J^dx-i\v6u\^ + J^dx 5-f\V {u + 6u)\^ + J^dx 5u ■ Vu (21) 

+ Yzi f ds\5j\^ + Y,Z(. [ ds5i~i = Q 
e=i ■^^'^i 1=1 -^^^e 
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In B, subtracting V • 7VU = from V • (7 + 5^)V{u + 6u) - gives the elliptic equation 

V ■[jV6u + 6-fV{u + 6u)]^0 iriB, (22) 



and then applying ( 18 ) for ifj - 5u yields 

/ dx7|V(5iip + / dx 5^\/ {u + 5u) ■ V 6u 

- I ds^5uV5u-'n+ I ds 5^ 5uV {u + 5u)^ ■ n 

= J ds 5u 6j 

where the second equality holds true by the definition of the perturbed normal component 
of boundary current density j + 5j = (7 + 6^)V{u + 5u) ■ n. Substituting back to (21 ) yields 

dx 57I Vn|^ + J^dx + 5^)\/6u ■ Vu + 

+ [ ^s6fj+ f ds&il6j^0, (23) 

while applying the perturbations to the electrode potential boundary condition Q gives 
5u = SUi - Zi6j, and therefore the last integral term becomes 

f dsduSj ^ [ i^Ui - ze5j)6j 

= f ds5j-Y.ze f ds\6j\^ 

where the last equation is due to the following lemma. 

Lemma 3.1 The perturbations in electrical admittivity 7 7 + 5'j, and induced electric 
potential in the interior of the domain u u + 6u give rise to a perturbation in the boundary 
current density with vanishing integral 



J^^ds6j{x) ^0, x€dB. 



Proof From the Neumann boundary condition ^ the current applied at the £'th electrode 
satisfies 

le - dsjS/u-n- / dsj. 
Keeping fixed before and after effecting the perturbations gives 

h- j ds {■^ + 5^)V {u + 5u) ■ n - i ds{j + 5j). 

Splitting the last integral, equating the right hand sides of the two equations above, and 
recalling from ([2]), that (j(x) + 6j{x)) = for X e 55 \ Te yields the result. 
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Effectively equation (p3|) reduces further to 

L 

^ da; 57|Vnp + / dx (7 + 57)V(5n • Vn + 



dx5^\Vu^+ f dx (7 + 57)V(5n • Vn + ^ Z£ f dsSjj^O, (24) 



and using once again the perturbed Robin condition the last integral simplifies further to 
Y^zg ds6jj = Y^SUi ds j - ds6u-j 



L _ r 

= Y h^Ui - ds6u-j 

L _ n _ 

= Yi h^Uf - / ds 7 5u • Vu • n 

£=1 -'^^i 



Now, consider the adjoint field problem (l7|)-([9|) subject to a current I™ = I*^. Then by the 
properties of the complete electrode admittance operator A^^z it is easy to show that the 
adjoint solution v{'j,I^) coincides with u{j,I'^). Applying the divergence theorem to the 
adjoint field equation Q gives 

y~ dxjVu-Vdu- J ds5vqVu-'n.- J dsSuj. 

Prom the above the perturbed power conservation law finalizes to 

Y h^Ui ^ - J^dx 57I Vnp - ^ dx djVSu • Vn - ^ dx (7 - 7) V5u • Vn. (25) 

Notice that for the purely real conductivity case, i.e. the cases of electrical resistance tomog- 
raphy where to - 0, the third term on the right hand side vanishes and the above collapses 



to the formula provided in 36 



Lemma 3.2 If the applied currents are purely real, the perturbed power conservation law 
(25) simplifies to 

Y h^Ui = - / dx 57 Vn • Vn - / dx 5jV5u ■ Vn. (26) 



Proof Consider applying the diverge theorem to (22) for a test function tp - u and to the 
adjoint pde ([I]) for = Su. Then upon subtracting the later from the former yields, 

Ylg^Ue = - J^dxSj\\/u\'^ - J^dxSj\/5u-\/u- J^dx{j -j)\/5u-\/u 

= / ds75nVn-n- / ds n(7V(5n + (57V(n + 5n)) • n 

= ds6uj- J~ dsu6j 

= / ds {6Ui- Zi6j)j - / ds {Ue- zij)6j, 
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where the last equation is due to lemma (3.1). Similarly, from the diverge theorem to (22) 
with f - u and to ([T]) with ip = 5u one obtains 



IB 



dx 5^\/u ■ \/u - / dx 5"f\/5u ■ \/u 



IB 



- J ds^duVu-Ti- J ds n(7V(5u + 57V(u + 5n)) • n 
= J ds6uj- J dsu5j 

= / ds {5Ui>- zi5i)j - / ds {Ue- Zij)6j ^Y,Ii6Ue. 

From the above the result follows in the case where = i.e. the imaginary component of 
the currents is zero. 

For simplicity we assume the case of real excitation currents. For a current pattern /, 
let 7p, 7 e L°°(B), the states of the model before and after the admittivity perturbation so 
that the change on the potential of the tth electrode is 

SUe{I)^Ui{^,I)-Ue{jp,I), 



and evaluate equation (25) for some pair drive current patterns that satisfy the constraint 
nbn. Let Hd, ^J'm ^ as in (12) some discrete patterns of zero sum, and define the currents 



1 = a/xrf. 



Suppose the currents are applied to the model with known admittivity 7p, and then to that 
compute the difference as 



of the unknown 7, giving rise to U{'^p,P) - A^^^^P , and U{'y,P) - A^^^P, from which we 



6UiP)^Uij,P)-Ui^p,P), 
for t - {d,m,c}. Based on the linearity of the admittance operator we deduce that 

6u{P) = A-^M" + n - K^ii' + n, 

6U{P) = A-^.P - A'l.P, and 5C/(/'") = A-^J"^ - A-lJ"^. Evaluating the left hand side 



of (25) for the three current patterns yields 

L 



e=i 1=1 e=i 

It is worth noticing that only SW^ are realistically measurable, since data acquisition occurs 
only under the direct patterns and borrowing the reciprocity result (11) for Iq = 1 gives 

j:p,6u^-j:if5u^-j:ir5ur = 2(5<, 



(27) 



£=1 



Expanding the corresponding right hand sides from ( 26 ) yields 

Y,P^5Ut-Y,lf5Uf-Y,lT5Ur^-2 f dx 6j Vu{PyVu{r')-2 f dx6-fV5u{P)-Vu{r'), 
1=1 i=i 1=1 -^^ -'^ 
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where we have used u{^,I'^) - 71(7,/'^) + u{j,I"^) for the interior fields. Let the fc'th mea- 
surement be = IJ-mU and note that u{"fp,I'^) - v, for v the adjoint fields solution of 
In effect, substituting and simplifying yields 

6Ck-- f dx57 Vn(7p,/'^)-V7J(7p,/")- f dx 5j VSu{I'^) ■ S/v{-fp, I""). (28) 

We are now ready to tabulate our main result in the form of the following theorem. 

Theorem 3.3 (The forward EIT transform) Consider the complete electrode model of ^ - 
1^ on a simply connected domain B, and suppose assumptions^hold. Suppose further that 
the applied currents are purely real and that boundary measurements C, e C" are observed. 
If u is the direct solution of this problem and v the pertinent adjoint vector satisfying 
then for any prior admittivity guess 7p e L°°{B) with direct solution £{'yp), the data change 
6Ck the kth element of the residual 5C - C - Sijp) satisfies 

5Ck-- f dx 5-/\/u{-f) ■ S/v{jp), (29) 

where ^7 = 7 - 7p is the residual vector between the target solution and the initial-prior guess. 

Proof The result follows immediately by substituting 6u - u{'^)-u{^p) for all direct currents 



I'^ to the integral equation (28), and holds true for all admissible bounded perturbations (57. 
This completes the proof. 



We would like to note that, in the Appendix we provide an alternative derivation of (29) 
suggested to us by an anonymous reviewer based on a weak formulation of the problem. 

3.2 Generalization to Poisson's equation with mixed boundary conditions 

Although the complete electrode model is now widely used for EIT, our new model formu- 



lation in (29) as well as the image reconstruction method to be described next are easily 
amenable to treat more simplistic electrode models. In particular, we now show that the 
above result holds true for a more general setting of impedance imaging involving the Poisson 
equation with Dirichlet and Neumann boundary conditions and point electrodes [2] , [24j . In 
geo-electrical application one usually encounters the model 

V-[7(x,^)Vn(x,w)] = /(x), xeS, (30) 

with boundary conditions of the form 

a(x)7(x,a;)Vu(x,a;) • n + /3(x)u(x,a;) = 0, x e 9-6. (31) 

where a and /3 are functions defined on dB and are not simultaneously zero to thoroughly 
impose the boundary conditions. To consider problems with different types of boundary 
conditions on different regions of dB, the functions a and (3 are allowed to be discontinuous. 
Figure [2] shows a common geophysical problem associated with the model in (30)-(31). In 



this problem r„ is the interface between the earth and air where a zero current condition 
(/? = 0) holds. In the remaining boundary Tm = dB\Tn, the values a and f3 are appropriately 
chosen to model an infinite half-space fSS^. When the sources of current are far from r^, 
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Figure 2: Geophysical application problem setting. The current sources f(x) are applied 
though the borehole electrodes yielding electrode potentials Ui. Tm is the model termination 
boundary and r„ is the upper surface of the model domain B. 
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a zero potential condition (a = 0) may be used as an approximation to the infinite half- 
space [42 j. 

The electric potential measurements are collected through point-wise electrodes, contact 
impedances of which are effectively zero. The measurement points are for i - 1,2, ... ,L 
and the measured potential at every point is 

Ui^ f dxu(x)5(x-Xi), (32) 

where 5{.) denotes the Dirac delta function. Consider a perturbation 7 7 + (^7 in the 
additivity causing the potential perturbation u ^ u + 6u. Introducing these into (30)-(31) 
gives 



a{j + 6'j)V{u + 5u) ■ n + I3{u + 5u) - 0, 



on B, 
on dB. 



Expanding (32) and (33) and using (30)-(31) to simplify the resulting terms yields 



V • (S'jVu) + V • i'yVSu) + V • {5^V5u) = 0, 
a{6^Vu - n + j\/6u - n + 6^\/6u - n) + f3Su - 0, 



on B, 
on dB. 



(33) 
(34) 



(35) 
(36) 



Based on (32) a perturbation in the measurement at x^ can be written as a volume integral 

6Ui^ f dx5n(x)5(x-X£) (37) 

To proceed with finding a closed form for the measurement perturbation 5Ui, it is useful to 
define v, as the solution to the adjoint system 



V-(7Vw£) = 5{yL-xg), 
d^Vvii ■ n + [ivi - 0, 

from which it is easily inferred that vl satisfies 

V-(7VW) = (^(x-xf), 

wyVUl ■ + (iUl - 0, 



xeB, 
xe dB, 



xeB, 
xe dB. 



(38) 
(39) 



(40) 
(41) 



terms of the adjoint field as 



Using (37) and (40) we conclude that the perturbation to the residuals can be written in 

dx 5u{x) V • (7VW)- 



6Ue 



IB 



(42) 



The remaining derivation requires extensive use of the following identity derived from Green's 
theorem [30] for vector function ^ and scalar function t/j 



/ dx*-V'0+ / dx^V-*= / dsip^-n. 
Jb Jb JdB 



We begin by taking il) - 5u and ^ = 7Vt'^ in (42) to obtain 



SUf - - / dx^S/Vi ■ \/Su + / dsj6u\/v£-n. 
Jb JdB 



(43) 



(44) 
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Next using ip - and ^ - "f\/5u in the first term on the right hand side of (44), we have 

6Ui^- J dxvgV ■ {^V6u) - J' ds^V£V6u-n+ J ds^fduVvl- n. (45) 



Prom (35), V-(7V(5u) = -\J ■ {5^\Ju) -\J ■ {5^\J 5u) which we use in the first term on the right 



hand side of (^45^ to arrive at 



5U^-- J dxvfV ■ i^5^V{u + 5u)^ - J ds^f V(V5u--n. + J ds^SuVv^- n.. (46) 



Appealing once more to (43) with ip - V£ and - 6^\/{u + 6u) in the first term of (46) gives 



6U£ - J dx5^\/Vf -Vu + J dx5^Vvn-V6u 

-J ds {jV(V5u ■ n + SjVfiS/u ■ n + S^UfV^u ■ n - ^SuVvl ■ n) . 



(47) 



We now show that the surface integral term in (47) is zero. For this purpose we multiply 

(48) 



both sides of (41) by 5u to arrive at 

aj5uVv£ ■ n + fiSuUl - 



Using (36) to replace the term fidu in (48) results in 



- a{^-~iV(V6u ■ n + SjVf S/u ■ n + djvi'VSu ■ n - 'jSu\/Vi • n) = 0, on dB. 



(49) 



The parenthesized expression in (49) is the same as the surface integrand in (47). We 



partition the boundary dB into Fq, where a ^ and dB \ Fq, where a - 0. Clearly (49) 



results the inside bracket expression to vanish on Fq; . On the remaining surface dB \ Fq, 
that a = 0, we certainly have /3 ^ since a and f3 may not be simultaneously zero and using 



this fact in (36) and (41) would result in 5u - and Vi - which again make the inside 



bracket term zero. Therefore the surface integral in (|47j) vanishes both on Fq, and dB \ Fq, 
and therefore 

6Ue,^ j^dx5^Vv~rVu+ J^dx5-iVvi-V5u, (50) 
and thus by substituting for 5u in the second term we arrive at the result of the theorem 



4 High-order regularized regression 

Within the d- dimensional sphere S^^^^, the electric potential field in the interior of the 
domain admits a Taylor expansion 

n(7) = n(7p) + 9^n(7p)57 + 0(||<572||) 

hence to first-order accuracy this can be approximated by 

u{'y) ^u(j) ^ u{jp) + d^u(jp)5-f. (51) 
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Introducing the right hand side of (51) in the integral equation (29) gives 



dx 57V(u(7p) + dyu{jp)6'y) ■ V-u(7p), 



IB 

/ dx 6-fVu{jp) -Vvijp) - / dx 6jV{d.yu{jp)6j) -Vvi-fp), 



(52) 



where the first, hnear term, involves the definition of the Frechet derivative of the forward 
mapping as in ( |16[ ) [2], [32], and the second nonlinear term the differential operator d^u{jp) 
that provides a measure on local sensitivity of the potential in the interior of the domain to 



perturbations in electrical properties. From (51), (52), it is trivial to deduce that the linear 



approximation of the forward operator £^ as in ( 15 ), as proposed by Calderon in |8], effectively 



imposes a zeroth-order Taylor approximation on the electric potential 14(7) - u{jp). In turn 
this enforces d^u and higher-order derivatives to vanish everywhere in B, thus elliminating 
the nonlinear terms in (28) and (52). Let the linear operator d-y£ = ■ L°°{B) C™, and 
nonlinear, quadratic in (57, /C : L°°{B) C" defined by 



J 57 



dx 57 Vn(7p) • Vv(7p), 



/C(57 i - / dx 57 V3^ti(7p)(57 • Vt'(7p) 
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(53) 
(54) 



then the inverse problem can be formulated in the context of regularized regression based 
on the nonlinear operator equation 



K - J 5^ + IC6j + rj. 



(55) 



4.1 Numerical approximation 

Usually the EIT problem is approached with a numerical approximation method like finite 
elements, where the governing equations are discretized on a finite dimensional model of the 
domain, say Bh{n,N) comprising n nodes connected in elements [37]. For simplicity in 
the notation we assume linear Lagrangian finite elements and consider element-wise linear 
and constant basis functions for the support of the electric potential u and conductivity 7 
respectively, 



i=l 
[N 



■■Bh 



N 

1=1 



Xi 



Bh 



(56) 



where {(/>i}"^^ and {xi}]ti the respective bases in B^- Following the discretization of the 
domain into a finite number of elements, the basis functions . . . , (pn) in the expansion 
of the potential are assumed to belong in a finite, n-dimensional subspace of H^{B). For 
clarity in the notation, we keep u and 7 as the vectors of coefficients relevant to the respective 
functions as from now on we deal exclusively the numerical approximation of the problem. 



On the discrete domain the weak form of the operator equation (55) is approximated by 



, m 



(57) 



where e C is the kth measurement, jk the A;th row of the Jacobian matrix J that is the 
discrete form of d^£{'yp), e (^NxN -g ^^iq kth coefficients (Hessian) matrix derived from 
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/C in (54), the noise in the kth measurement and 6^ e C''^ the required perturbation 
in the admittivity coefficients. Let the additive noise be uncorrelated zero-mean Gaussian 
with diagonal covariance matrix Cn, with positive diagonal element then the data misfit 
function 

m 

Q{57)-E^k\SCk-jkSl-S7'K'6jf (58) 
fc=i 

can be used to define the regularized quadratic regression problem 



1 



57* = arg min £((57), ^{S'y) ^ -{Q{6-f) + ag{6-f)} 



(59) 



with g ■ M a convex differentiable regularization term. On the other hand, choosing 

to neglect the matrices K'^ yields the conventional misfit function 

m 

M5^)-j:^k\SCk-jkSl) , (60) 
fc=i 

often used in the context of regularized linear regression formulations. As shown in [l] the 
Jacobian matrix can be computed directly from (53) and (56) using numerical integration 
as 



'k,j - 



Zesupp(_Bj) I esupp(Bj) 



k^l,...,m, j = 1,...,N 



(61) 



with V the coefficients of the adjoint field solution corresponding to the A;th measurement, 
and supp(-Bj) the support of the jth element. To derive the respective element of K'^ we 
follow an approach similar to that of Kaipio et al. in [21] that is based on the Galerkin 
formulation of the problem. For this we choose {(pi, . . . , as a test basis for the potentials 
and by substituting into the variational form of the model we arrive at 

n n p L ^ \ ^ r 

dxjV4>i-y(t>j + Y.ze ds4>i(f)j]ui-Y,^e ds<piUe^O. 

i=l j=l -'^ 1=1 ~'^<^i 1=1 ~'^<=e 

Imposing the Neumann conditions for the applied boundary currents yields the additional 
equations 

h^-zeY,i ds(l)ijui + ze\Tei\Ue, £ = 1,...,L, 

with iFe^l the area of the £th electrode. In matrix form the electric potential expansion 
coefficients u e C" and the electrode potentials U e can be computed by solving the 
(n + -L) X (n + L) matrix equation 



'-12 



A'i2 A22 





u 









U 




/ 



(62) 



where 



An 



/ dxjVcpi ■ S/4>j + Y,ze di 



i,j = l,...,n 



A12 j,^ 

A22^,^ 



Zel^eel, 



ds < 



z = l,...,n, £=1,...,L, 



^1,...,L. 
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For a conductivity 7 and applied current /, let 



the solution of (62). Using 



the matrix differentiation formula, the partial derivatives with respect to the qth admittivity 
element are 



-A-i(7)9,,{A(7)}A-H7)/ 
-A-i(7)9,,{A(7)} 



where 



a^,{A(7)} = 5^jAii(7)}= r dxV0»-V</.j, q^l,...,N, 

J Bq 

as only the block An depends on admittivity. Separating the above as 



and evaluating the upper part for all elements in the model yields the required matrix in 
vector concatenation form 



d-yu{-i) = [ 5^1^(7) I dy,u{-i) I ... I d^^u{-i) ], 



(63) 



while d^^U are the elements of the Jacobian matrix J. Effectively the element of K*^ matrix 
is given by 

K^j^-J dxtl^j Vcpidy^ui ^ vi\/(t>i, k^l,...,m, r,j ^1,...,N 

J ksupp(Br) iesupp(Bj) 

with V the adjoint field corresponding to the kth measurement and d^^u the derivative of 
the kth direct field with respect to 7^. 

4.2 Newton's minimization method 



We propose solving the regularized problem (59) using Gauss-Newton's minimization method 
[18|. At a feasible point ^7^ the minimization cost function ^ is approximated by a second- 



order Taylor series 16 



e(<57) - ^i^lp) + ds.,^{Sip){Sj - Sjp) + -(57 - SjpYds^s^^SjpXS^ - 5jp), (64) 
where applying first-order optimality conditions ds-y£,{Sj) - yields the linear system 



Prom ( 58 ) , let the feth residual function be 

N 



N N 



j=i j=i 1=1 
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such that Q{6^) - \\r{6j)\\'^, then the cost gradient ds-y^i^p) and Hessian ds^Sj^i^lp) are 
expressed as 

dsjS-yii^lp) = dsyr{6-fpyds^r{6'yp) + aC:^^ 

for r((57) = [?"i((57), . . . , rm{S^)]' , and assuming a Tikhonov-type regularization function 
^((^7) = a6^'C^^6j, with C^^ positive semidefinite and a a positive regularization parame- 
ter. The Jacobian of the residual dsjr{djp) e C"^'*^ is then formed using the vectors 



N 



evaluated at 6^p like 

ds^r{5jp) ^[ds-yri{5-fp) \ ds^r2{5-ip) \ ... \ ds^rm{5lp)\ . 

If ds-yS-yCi^lp) is full rank and positive definite the solution can be computed iteratively using 
Newton's algorithm 

hp+i^ Sjp-ds^s^ii^^p)95ji{Sjp), p = 0,l,2,... (65) 

Using standard arguments from the convergence analysis of Newton's method on convex 
minimization it is easy to show convergence as in ^18^ , ^27j 

l(57p)>e(<57p+i)>lhll, \\Sl*-S7p\\>\\Sl*-6jp^i\\, p^O,l,..., (66) 

however a convergence in the sense of the discrepancy principle is more appropriate as the 
data are likely to contain noise 122]. 



Corollary 4.1 Initializing the quadratic regression iteration (65) with - yields a first 
iteration that coincides with the linear regularized regression estimator 

571 = (J'C^^J + aC-^y^3'C-^^5(: (67) 

Proof The proof is by substitution of the residual and its Jacobian at - into the 
expressions for the gradient and Hessian of the cost function. In particular for r(57o) = 

—1/2 —1/2 

5C, and ds'yr{5^Q) - J, iteration (65) yields the result. 



Combining the convergence remarks of (66) with the corollary above, we assert that for 
p> 1 the quadratic regression iterations should converge in a solution whose error does not 
exceed that of the linear regression problem (17). Suppose now that at a certain iteration p 
the value of the residual r{6^p) converges to the level of noise ||r/||. Then according to the 
discrepancy principle one updates the admittivity estimate as 7^+1 = 7p + ^7^ and thereafter 
the definitions of J and K'^, and then proceeds to the next iteration. Effectively, the resulting 
scheme can be expressed as a Newton-type algorithm. 
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1. Given data C e with noise level ||?7|| and a finite domain with unknown admit- 
tivity 7* e 

2. Set q - 0, choose initial admittivity distribution 70, 

3. For q - 1,2, . . . (Exterior iterations) 

4. Compute data 6C = C -<?(7g-i), and matrices J e C""'^^, K'^ e C^'^^, for 7g_i, and 
/c = 1, . . . , m, 

(a) Set p - 0, 5^p - 0, 

(b) For p - 1,2, . . . (Interior iterations) 

(c) Compute update 

S-/p = 6jp-i - Tpds^s^i{5jp-i) ds^i{S-fp-i), Tp > 0, 

(d) End p iterations 

(e) Compute update 

7g = 7g-l + Tq^lp, Tq > 0, 

5. End q iterations 



In performing the outer iterations, a complication will likely arise in that a certain update 
admittivity change 5'~fp may cause the real and/or imaginary components of 7^+1 to become 
zero or negative. This of course violates a physical restriction on the electrical properties of 



the media, and the solution cannot be admitted. For this reason the problem of (59) should 
be posed as a linearly constrained problem 

(57* = arg min ^(57), 

7,><57 

at each 7^. A convenient heuristic to prevent this complication is by adjusting the step sizes 



Tq, Tp until the above inequality is satisfied 37 , 43 . Note also, that the above methodology 



makes no explicit assumptions on the type of the regularization functional ^(7), aside its 
differentiability, thus we anticipate it can be also be implemented in conjunction with total 
variation and £i-type regularization as well as the level sets method [llj . 

5 Numerical results 

To test the performance of the proposed algorithm we perform some numerical simulations 
using two-dimensional models, although the extension to three dimensions follows in a trivial 
way. In this context we consider a rectangular domain B - [-16,16] x [0,-32] c M^, with 
L = 30 point electrodes attached at its boundary in a borehole and surface arrangement as 
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shown in figures |6] and [7j As a first test case the domain is assumed to have an unknown target 
conductivity 7* whose real and imaginary components are functions with respective bounds 
1.46 < a* < 5.60 and 0.74 < uje* < 3.90. To compute the measurements we consider 15 pair 
drive current patterns I'^, d - 1, . . . ,L/2, yielding a vector of m = 390 linearly independent 
voltage measurements ^ ^ C"^. The forward problem is approximated using the finite element 
method outlined in the previous section, and to the measurements we add a Gaussian noise 
signal of zero mean and positive definite covariance matrix = 10"^ max |C| I, where I is the 
identify matrix. For the forward problem we use a finite dimensional model Bf comprising 
n - 1701 nodes connected in = 3144 linear triangular elements. All other computations are 
performed on a coarser grid Bi with n = 564 nodes and A^ = 1038 elements. The two finite 
models are nested, hence for any function 7 approximated on Bi with expansion coefficients 
7j there exists a projection 7^ = Hji, mapping it onto Bj. To reconstruct the synthetic data 
we assume an initial homogeneous admittivity model 70 = 3.90 + 2.40* which coincides with 
the mean value of 7*, a methodology adopted from [19| . 

At the initial admittivity guess 70 we approximate the potential u(^*) using the zeroth- 
order and first-order Taylor series ^^(70) and ^(70) + d^u(5jQ)(^ - 70) respectively. The 
normalized approximation errors are illustrated at the top of figure |3] next to those of the 
error in the induced potential gradient as this is involved in the computation of the K'^ 
matrices for k - 1, . . . , 390. The results show that the linear approximation sustains a smaller 
error in both quantities and at all applied current patterns. In the same figure we also plot 
the measurement perturbations 5^ - ^-^(70) versus the linear and the quadratic predictions 
to demonstrate that the proposed quadratic regression will fit the noisy measurements at a 
smaller error. In particular, the quadratic and linear misfit cost functions in (58) and (60) 
are evaluated at Q{6^o) - 0.06 and A{Sjq) - 0.13, where ^70 = 7* -70- Notice the impact of 
the second-order term, that brings the norm of the data misfit to about half of that of the 
linear case. 

To reconstruct the admittivity function we implement the proposed iteration (65) using 
a precision matrix C^^ = R'R, where R e M.^^^ is a smoothness enforcing operator. In the 
numerical experiments we use two different values of the regularization parameter in order 
to investigate the performance of the scheme at different levels of regularization. Using 
a - 5 X and a = 5 x 10~^, we execute two exterior GN iterations each one incorporating 
two inner Newton iterations after which the algorithm converged to an error value just above 
the noise level. The error reduction is illustrated by the graphs of figures [4] and [5| showing 
a significant reduction in both the misfit error Q{6^p) and the image error \\S'y* - Il6'yp\\ 
respectively for = 0, 1, 2 for the first and second GN iterations, i.e. q - 1,2. Each exterior 
iteration was initialized with ^70 = hence we can regard 5ji as the Tikhonov solution (67) 
and 572 the quadratic regressor after two iterations ( 65 ) . To aid convergence a backtracking 
line search algorithm was used where the optimal step sizes for each iteration, interior and 
exterior. The computational time required to assemble the Jacobians J e ^390x1038 
about 0.34 s, while each of the 390 matrices K'' e ^1038x1038 ^q^j^ about 4.75 s and then 
each iteration about 12 s depending on the line search. These times are based on running 
Matlab fsT] on a machine with a dual core processor at 2.53 GHz. Despite the substantial 
computational overhead, the method can be appealing in the cases where the inverse problem 
is heavily underdetermined with only a few measurements. Moreover, the assembling of the 
K matrices is well suited for parallel processing. The images of the reconstructed admittivity 
perturbation at each iteration are plotted in figures [6] (real component) and [T] (imaginary 
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Figure 3: At the top row, the normahzed errors in the electric potential field approximation 
and its gradient, assuming zeroth-order (dashed line with + markers) and first-order (solid 
line with x markers) Taylor series approximations of u{j*,I'^) direct fields. In both cases 
the errors with the linear approximation are lower. Second row, the quality of the linear 
and quadratic approximations in predicting the nonlinear change in the boundary data S(^. 
The solid line denotes S^, the dashed j'^dj and the dotted j'^Sj + Y.ili S^K^S^, over the 
interval i = 150, . . . , 220. The corresponding data misfit norms are 0.057 for the quadratic 
approximation Q{Sj*) and 0.123 for the linear A(d'j*), assuming no additive noise. With 
the prescribed additive noise these values change to 0.062 and 0.126 respectively. 
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Figure 4: Indicative convergence of the proposed method, in terms of minimizing the 

quadratic misfit error Q{5^p) for two different values of the regularization parameter a. 
Left the results during the first exterior iteration q = 1, and right the corresponding values 
for q = 2. In these results, ^70 = 0, (^71 coincides with the Tikhonov solution, and (^72 is 
the regularized quadratic regression solution. Notice that the quadratic regression solution 
has lower data misfit errors in both GN iterations. Between the first and second exterior 
iteration the admittivity increment was scaled to preserve positivity, hence the apparent 
discontinuity in the error reduction. 




Iteration index p 



Iteration index p 



Figure 5: Indicative convergence of the proposed method, in terms of minimizing the the 
image error Wdj* - TlS'fpW for two different values of the regularization parameter a. Left the 
results during the first exterior iteration q= 1, and right the corresponding values for q = 2. 
In these figures ^70 = 0, Sji coincides with the Tikhonov solution, and 6^2 is the regularized 
quadratic regression solution. Notice that the quadratic regression solution maintains lower 
image errors at each external iteration. Between the first and second exterior iteration the 
admittivity increment was scaled to preserve positivity, hence the apparent discontinuity in 
the error reduction. 
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Figure 6: At the top, the simulated target conductivity a* profile on Bf, as used in the first 
test example and the arrangement of the electrodes. In the second row, from left to right, 
the respective images resulted from first exterior iteration using a - 5 x 10"^, namely the 
real components of 70 + t5ji, 70 + tSj2, and 71 on Bi. Similarly at the bottom row, the 
respective images from the second exterior GN iteration, 71 + tS'Ji, 71 + r(572, and 72, using 
the same value of a. 
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Figure 7: At the top, the simulated target scaled permittivity uje* profile on Bf, as used in 
the first test example and the arrangement of the electrodes. In the second row, from left to 
right, the respective images resulted from first exterior iteration using a - 5 x 10"^, namely 
the imaginary components of 70 + t5^i, 70 + t572, and 71 on Bi. Similarly at the bottom 
row, the respective images from the second exterior GN iteration, 71 + t6^i, 71 + t5j2, and 
72, using the same value of a. 
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component) below their respective target images for comparison. As the error graphs clearly 
indicate, the reconstructed images show a profound quantitative improvement in spatial 
resolution, with the regularized quadratic regression solution 5^2 to outperform the Tikhonov 
solution in both Gauss-Newton iterations. Notice however, that in the exterior iteration 
we scaled the increment 5^p by Tq in order to preserve the positivity new admittivity estimate. 
This scaling, if < 1 tends to increase the data misfit errors, hence one can observe some 
discontinuities in the error reduction from q-lioq-2'm. the plots of figure |4j Similarly for 
the graphs of the image error in figure [Sj although this time the correction works out to the 
improvement of the errors as the target images are by definition positive. For completeness, 
the step sizes used in the image reconstructions of figures [6] and [7] are Tp=i = 1, rp=2 = 0.3, 
and Tq^i - 0.78 for the first cycle of iterations and rp=i = 1, rp=2 = 0.38, and rg=i = 1 for the 
second. 

As a second example we consider a purely conductive case, i.e. w = 0, aiming to re- 



construct the target conductivity function appearing at the top of figure 10 Once again 



synthetic data are simulated, using the same current and measurement patterns as in the 
previous case. After computing the measurements C and introducing some zero mean Gaus- 
sian noise using the noise covariance covariance matrix = 10"^max|(^|I we formulate the 
inverse problem at a homogeneous background conductivity do, the best homogeneous fit 
of the data, regularization matrix R, and a parameters equal to 10^^ and 10~^. To aid 
comparison with the convergent results for the complex admittivity case we implement the 
algorithm for two interior and two exterior iterations, for each of the regularized problems. 
The graphs of the data misfit and image errors are illustrated in figures [8] and [9} The graphs 
show a convergence pattern similar to that of the complex case, for both values of the regu- 
larization parameter. Also at the initial reference point the linear and quadratic data misfit 
functions obtain values A((^cJo) = 0.24 and Q{5aQ) - 0.09, demonstrating once again that the 
contribution of the quadratic term can be significant if the reference point is not sufficiently 
close to the solution. In terms of its computational cost, implementing the algorithm for the 
purely real admittivity has brought the processing time to about a half of that consumed 
for the complex case. The reconstructed images presented in figure 10 correspond to the 
various conductivity updates as computed for two exterior and two interior iterations, with 
a - 10^^. Initializing with 5aQ - the second row, from left to right, shows the conductivity 
updates after each interior iteration for the first exterior GN iteration, and the bottom row 
the respective images from the second exterior iteration. The step sizes used in these results, 
as computed by the line search algorithms are rp=i = 1, Tp=2 = 0.3, and rq=i = 0.62 for the 
first cycle of iterations and rp=i = 0.62, rp=2 = 0.24, and rq=i = 1 for the second. Moreover 
at the beginning of the second GN iteration the misfit functions have been computed at 
K{5ai) = 0.034 and Q{5ai) = 0.007. 

When the noise level in the data is approximately known, solving the nonlinear EIT 
problem one typically performs a number of GN iterations until convergence is reached 
in the sense of the discrepancy principle [27] , [43] . In our results we implement only two 
exterior GN iterations, i.e. q= 1,2, each one encompassing two interior iterations, in order 
to demonstrate the observed reduction in the image and data misfit errors. Consequently, 
by virtue of the convergence properties of the Newton algorithm, it is straightforward to 
state that the quadratic regression solution will sustain a smaller error for any number of 
GN iterations [3], and will thus converge to the solution faster. On the other hand, a serious 
bottleneck of the second, respectively higher-order, formulation is the computational demand 
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Figure 8: Indicative convergence of the proposed method, in terms of minimizing the 
quadratic misfit error Q{6^p) for two different values of the regularization parameter a. 
Results are from the second test case with simulations at dc conditions cj = 0, hence the 
admittivity is purely real, i.e. 7 = 0". Left the results during the first exterior iteration q = l, 
and right the corresponding values for q-2. In these results, ^70 = 0, S^yx coincides with the 
Tikhonov solution, and ^72 is the regularized quadratic regression solution. Notice that the 
quadratic regression solution has lower data misfit errors in both GN iterations. Between the 
first and second exterior iteration the admittivity increment was scaled to preserve positivity, 
hence the apparent discontinuity in the error reduction. 



to compute the K matrices. In this sense the method is more suited to the cases where high 
performance computing is available, or when the number of data m is fairly small. 

6 Conclusions 

This paper proposes a new approach for the inverse impedance tomography problem. Based 
on a a power perturbation approach we derive a nonlinear integral transform relating changes 
in electrical admittivity to those observed in the respective boundary measurements. This 
transform was then modified by assuming that the electric potential in the interior of a 
domain with unknown electrical properties can be approximated by a first-order Taylor 
expansion centered at an a priori admittivity estimate. This framework yields a quadratic 
regression problem which we then regularized in the usual Tikhonov fashion. Implementing 
Gauss-Newton's iterative algorithm we demonstrate that the method quickly converges to 
results that outperform those typically computed by applying the algorithm on the linearized 
inverse problem. An important shortcoming of this approach is the computational cost of 
computing the second and higher derivatives, as they require the assembly of large dense 
matrices of dimension equal to that of the parameter space. A possible remedy to this can 
be found in model reduction methods |28j. Another interesting extension is to consider a 
reformulation of the inverse problem in terms of some surrogate parameter functions, e.g. 
the logarithm of the admittivity, in a way that preserves the necessary positivity on the 
electrical parameters. 
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Figure 9: Indicative convergence of the proposed method, in terms of minimizing the the 
image error ||57*-n(57p|| at each interior iteration for two different values of the regularization 
parameter a. Results are from the second test case with simulations at dc conditions a; = 0, 
hence the admittivity is purely real, i.e. 7 = <t. Left the results during the first exterior 
iteration g = 1, and right the corresponding values for q = 2. In the figures ^70 is the initial 
homogeneous guess, (^71 coincides with the Tikhonov solution, and ^72 is the regularized 
quadratic regression solution. Between the first and second exterior iteration the admittivity 
increment was scaled to preserve positivity, hence the apparent discontinuity in the error 
reduction. 



Appendix 



Here we present an alternative approach to the derivation of (29) suggested to us by one 



of the anonymous reviewers. Rather than relying on the conservation laws as was the case 
for our approach, the one presented below is based more on the use of variational methods 
applied to both the forward and adjoint problems. From the weak form ( 19 ) for ip - v{"fp, J™") 
and ^ - Vf(7p,/™-) assuming u{^,I'^), Ug{^,I'^) then 

[ dx7Vu(7,/'^)-Vt;(7p,/'") 

Repeating for a model with u{^p,I'^), Ue{^p,I'^) yields 
f dx7pVn(7p,/°')-Vr;(7p,/"') 

^t^i' L d.(u(7„/'^)-C/K7p,/'))(^(7p,0-m,0)-E^'^, 
and thus by subtracting and inserting ± J^dx"fp\/u{'j,I'^) ■ V'u(7p,/™') one arrives at 
0= f dx{^-jp)Vu{j,I'')-Vvhp,in+ f dxjpV{u{j,I'')-u{jp,I''))-v{jp,I^) 

J B J B 

+ i ds{u{^p,i')-u,{^p,i'')){v{^p,n-v,{jp,n)- 
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Figure 10: Simulated and reconstructed admittivity functions, for the second test case at 
direct current conditions. Top row, the target simulated conductivity a* discretized in Bf. 
Below from first exterior iteration using a = 10"^, namely ao + T5ai , ao + T6a2, and ai on Bi. 
Similarly at the bottom row, the respective images from the second exterior GN iteration, 
fji + T5ai, ai + Tda2, and cr2, using the same value of a. 
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Similarly from the weak form of the adjoint problem assuming ^(7^,/™) and Vi^jp, I"^) for 
ip = ti(7,/'^) and = Ueij,!'^) and ip = u(jp,I'^) and = Ui{jp,I''-) we get 

|:([/(7,/'') - t/(7p,/''))iF= rdx7pV(u(7,/^)-^x(7p,/'^)).Vt;(7p,/^) 
£=1 -^^ 



thus combining the last two relations yields the result (29). 
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